#include <stdio.h>
#include <stdlib.h>
#include <time.h>

/*	LU_decomposition.cpp
*	This program is part of the °Simple-Benchs° Benchmark
*	The complexity of this program is O(n³) where n is the Matrix's size
*
*	
*	To run the program:	./program MATRIX_SIZE MAX_VALUE 
*/


/* Array initialization. Complexity O(n²) where n is the matrix's sizeof*/
void initArray(float **vet, int size, int max_num){
  int rows, columns;
  for (rows = 0; rows < size; rows++){
    vet[rows] = (float*)malloc(size*sizeof(float));
    for (columns = 0; columns < size; columns++){
      vet[rows][columns] = (float)(rand() % max_num + size);
    }
  }
}

/* Program standart output. Display a Matrix.  Complexity O(n²) where n is the matrix's size */
void printArray(float **vet, int size){
  int rows, columns;
  printf("\n\t============ Matrix ============\n");
  for (rows = 0; rows < size; rows++){
    for (columns =0; columns < size; columns++){
      printf("\t%f", vet[rows][columns]);
    }
    printf("\n");
  }
  printf("\n");
}

/*	Initial check to verify whether the amount of arguments is correct or not.  Complexity O(1)	*/
short Initial_check(int inputs){
  if (inputs != 3){
    fprintf(stderr,"Wrong number of argments. Syntax: ./program (Matrix's size) (Max element)\n");
    return 1;
  }
  return 0;
}

/*	Auxiliar function to free all blocks in heap. Complexity O(n) where n is the amount of rows in the square matrix */
void free_all(float **matA, float **matB, float **matC, int size){
  int rows;
  for (rows = 0; rows < size; rows++){
    free(matA[rows]);
    free(matB[rows]);
    free(matC[rows]);
  }
  free(matA);
  free(matB);
  free(matC);
}


/*	Perform the LU decomposition using the Crout algorithm. Complexity O(n³) defined by the nested loops starting the line 71, 
where n is the number of rows in the input square matrix */
void luSolver(float **A, float **L, float **U, int n){
	int rows, columns, count;
	float sum = 0;
 
	for (rows = 0; rows < n; rows++){
		U[rows] = (float*)calloc(n,sizeof(float));
		L[rows] = (float*)calloc(n,sizeof(float));
		U[rows][rows] = 1;
	} 
	for (columns= 0; columns < n; columns++){
		for (rows=columns; rows< n; rows++){
			sum = 0;
			for (count = 0; count < columns; count++)
				sum = sum + L[rows][count] * U[count][columns]; L[rows][columns] = A[rows][columns] - sum;
		}
		for (rows=columns; rows< n; rows++){
			sum = 0;
			for(count = 0; count <columns; count++)
				sum = sum + L[columns][count] * U[count][rows];
			if (L[columns][columns] == 0){
				printf("Opssss .... L-Matrix's determinant is getting close to 0. It is not possible to divide by 0 yet.. So I going to kill myself ! ! !\n\n\n");
      		free_all(A,L,U,n);
				exit(1);
			}
			U[columns][rows] = (A[columns][rows] - sum) / L[columns][columns];
		}
	}
}

int main(int argc, char** argv){

	/* Checking arguments */
	if (Initial_check(argc)) return 1;

	/* Retrieve problem size. */
	int size = atoi(argv[1]);
  int max_num = atoi(argv[2]);

 	/* Seed initialization. */
 	srand((unsigned)time(NULL));

 	/* Variable declaration/allocation/initialization. */
 	float **matrix = (float**)malloc(size*sizeof(float*));
 	float **L = (float**)calloc(size,sizeof(float*));
 	float **U = (float**)calloc(size,sizeof(float*));
 	initArray(matrix, size, max_num);
  	
 	/* Display matrix A from A = LU */  	
 	printArray(matrix, size);
  	
 	/* Running the LU decomposition */
 	luSolver(matrix,L,U, size);

 	/* Display matrix U and L from A = LU */  	
 	printArray(L, size);
 	printArray(U, size);

 	/* Be clean. */
 	free_all(matrix,L,U,size);
  return 0;

}